clear

import delimited "U:\school testing\data_all.csv"

gen date1= date(date, "YMD")

egen pannel_var= group(related_organization poc_test_type reason_for_poc_test test_type account_role__c test_result)

gen cases_norm= (totalconfirmedcases/ county_population)
encode school_type ,gen(school_category)
xtset pannel_var date1

**Jan to feb28
xtreg positivity_rate cases_norm consent_rate_cumsum black hispanic white other i.rural i.title1 i.charter i.school_category i.is_low if school_level== "elementary" & date1 >= date("2022-01-01", "YMD") & date1 <= date("2022-02-28", "YMD")

xtreg positivity_rate cases_norm consent_rate_cumsum black hispanic white other i.rural i.title1 i.charter i.school_category i.is_low if school_level== "high" & date1 >= date("2022-01-01", "YMD") & date1 <= date("2022-02-28", "YMD")

**Mar1 to end
xtreg positivity_rate cases_norm consent_rate_cumsum black hispanic white other i.rural i.title1 i.charter i.school_category i.is_low if school_level== "elementary" & date1 >= date("2022-03-01", "YMD")

xtreg positivity_rate cases_norm consent_rate_cumsum black hispanic white other i.rural i.title1 i.charter i.school_category i.is_low if school_level== "high" & date1 >= date("2022-03-01", "YMD")

************************************************************

**Propensity score matching

clear

import delimited "U:\school testing\Revision-Dec 2023\Jan_2024\propensity_match_data.csv"

gen date1= date(date, "YMD")

gen surv_date= date(surv_policy_go_live_date, "YMD")

gen weeks_from_surv_start= round((date1 - surv_date)/14)

egen pannel_var= group(related_organization poc_test_type reason_for_poc_test test_type account_role__c test_result)

gen cases_norm= (totalconfirmedcases/ county_population)

encode school_type ,gen(school_category)

xtset pannel_var weeks_from_surv_start

teffects psmatch (positivity_rate) (group full_elementary_porportion full_middle_porportion full_high_porportion full_district_porportion consent_rate_cumsum black white hispanic other rural title1 charter masking_policy school_category), gen(match) nneighbor(1)

predict ps0 ps1, ps // The propensity scores for each group.
predict y0 y1, po // The potential outcome estimated for each group.
predict te // The treatment effect estimated for each group.

************************************************************
**GEE Regression

clear
import delimited "U:\school testing\Revision-Dec 2023\Jan_2024\target_trial_data_jan_19_1.csv"

gen date1= date(date, "YMD")

gen tts_date= date(tts_policy_go_live_date, "YMD")
gen surv_date= date(surv_policy_go_live_date, "YMD")

gen weeks_from_tts_start= round((date1 - tts_date)/14)
gen weeks_from_surv_start= round((date1 - surv_date)/14)

egen pannel_var= group(related_organization poc_test_type reason_for_poc_test test_type account_role__c test_result)

gen cases_norm= (totalconfirmedcases/ county_population) * 100000

encode school_type ,gen(school_category)

xtset pannel_var weeks_from_surv_start

gen weeks_from_surv_start_pos = weeks_from_surv_start + 13


*Elementary - surv_init

xtgee positivity_rate i.group cases_norm full_elementary_porportion i.weeks_from_surv_start_pos if school_level== "elementary" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >=0), corr(exchangeable) robust i(pannel_var)

*middle - surv_init

xtgee positivity_rate i.group cases_norm full_middle_porportion i.weeks_from_surv_start_pos if school_level== "middle" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >= 0), corr(exchangeable) robust i(pannel_var)

*high - surv_init

xtgee positivity_rate i.group cases_norm full_high_porportion i.weeks_from_surv_start_pos if school_level== "high" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >= 0), corr(exchangeable) robust i(pannel_var)

*combine - surv_init
keep if inlist(school_level, "all", "elementary/middle", "middle/high")

xtgee positivity_rate i.group cases_norm full_district_porportion i.weeks_from_surv_start_pos if (weeks_from_surv_start <= 6 & weeks_from_surv_start >=0), corr(exchangeable) robust i(pannel_var)

******************************************************************************
**sensitivity analysis - excluding 242 switchers


drop if (surv_date > date("2021-09-19", "YMD")) & (surv_date <= date("2021-12-12", "YMD"))

*elementary
xtgee positivity_rate i.group cases_norm full_elementary_porportion i.weeks_from_surv_start_pos if school_level== "elementary" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >=0), corr(exchangeable) robust i(pannel_var)

*middle

xtgee positivity_rate i.group cases_norm full_middle_porportion i.weeks_from_surv_start_pos if school_level== "middle" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >= 0), corr(exchangeable) robust i(pannel_var)

*high

xtgee positivity_rate i.group cases_norm full_high_porportion i.weeks_from_surv_start_pos if school_level== "high" & (weeks_from_surv_start <= 6 & weeks_from_surv_start >= 0), corr(exchangeable) robust i(pannel_var)

*combine

keep if inlist(school_level, "all", "elementary/middle", "middle/high")

xtgee positivity_rate i.group cases_norm full_district_porportion i.weeks_from_surv_start_pos if (weeks_from_surv_start <= 6 & weeks_from_surv_start >=0), corr(exchangeable) robust i(pannel_var)



******************************************
**TTS

clear 

import delimited "U:\school testing\data_all.csv"

drop if reason_for_poc_test != "Test to Stay"

gen tts_date= date(tts_policy_go_live_date, "YMD")
gen surv_date= date(surv_policy_go_live_date, "YMD")

gen weeks_from_tts_start= round((date1 - tts_date)/14)
gen weeks_from_surv_start= round((date1 - surv_date)/14)

egen pannel_var= group(related_organization poc_test_type reason_for_poc_test test_type account_role__c test_result)

gen cases_norm= (totalconfirmedcases/ county_population) * 100000

encode school_type ,gen(school_category)

xtset pannel_var weeks_from_tts_start

gen weeks_from_tts_start_pos = weeks_from_tts_start + 5


*elementary
xtgee surv_positivity_rate cases_norm full_elementary_porportion consent_rate_cumsum black white hispanic other i.rural i.title1 i.charter i.school_category i.weeks_from_tts_start_pos if school_level== "elementary" & (weeks_from_tts_start <= 6 & weeks_from_tts_start >=-2), corr(exchangeable) robust i(pannel_var)


*middle
xtgee surv_positivity_rate cases_norm full_middle_porportion consent_rate_cumsum black white hispanic other i.rural i.title1 i.charter i.school_category i.weeks_from_tts_start_pos  if school_level== "middle" & (weeks_from_tts_start <= 6 & weeks_from_tts_start >= -2), corr(exchangeable) robust i(pannel_var)


*high
xtgee surv_positivity_rate cases_norm full_high_porportion consent_rate_cumsum black white hispanic other i.rural i.title1 i.charter i.school_category i.weeks_from_tts_start_pos if school_level== "high" & (weeks_from_tts_start <= 6 & weeks_from_tts_start >= -2), corr(exchangeable) robust i(pannel_var)


*combine
keep if inlist(school_level, "all", "elementary/middle", "middle/high")

xtgee surv_positivity_rate cases_norm full_district_porportion consent_rate_cumsum black white hispanic other i.rural i.title1 i.charter i.school_category i.weeks_from_tts_start_pos if (weeks_from_tts_start <= 6 & weeks_from_tts_start >=-2), corr(exchangeable) robust i(pannel_var)



**Plotting TTS -

clear
import excel "U:\school testing\Revision-Dec 2023\Jan_2024\resultst_jan_19_2024\GEE_TTS.xlsx", sheet(plot) firstrow case(lower)


gen time_pp= time + 0.5
gen time_p= time + 0.2
gen time_m= time - 0.2
gen time_mm= time - 0.5

twoway rcap lb1 ub1 time_mm, lcolor(black) || scatter point1 time_mm, mcolor(pink) || rcap lb2 ub2 time_m, lcolor(black)|| scatter point2 time_m, mcolor(ebblue)|| rcap lb3 ub3 time_p, lcolor(black)|| scatter point3 time_p, mcolor(lavender) || rcap lb4 ub4 time_pp, lcolor(black)|| scatter point4 time_pp, mcolor(brown) yline(0, lcolor(maroon) lpattern(dash)) xline(0, lcolor(maroon) lpattern(shortdash)) title("Bi-weekly In-School Surveillance Testing Positivity Rate Stratified by School Type After Initiating Test to Stay", size(small)) xlabel(-4(2)12) xtitle("Weeks") ytitle("Positivity Rate") ysize(60) xsize(85) scale(0.9) legend(pos(6) col(1) order(2 4 6 8 10) label(2 "Elementary School") label(4 "Middle School") label(6 "High School") label(8 "Multiple Program School"))


